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Abstract. The hull-gradient method is used to determine the critical 
threshold for bond percolation on the two-dimensional Kagome lattice 
(and its dual, the dice lattice). For this system, the hull walk is 
represented as a self-avoiding trail, or mirror-model trajectory, on the 
(3,4,6,4)-Archimedean tiling lattice. The result pc = 0.524405 3 ± 
0.000 000 3 (one standard deviation of error) is not consistent with 
previously conjectured values. 
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1. Introduction 

The Kagome lattice (Fig. 1) is one of the fundamental lattices of two-dimensional 
percolation, as well as many other two-dimensional lattice problems. It is one of 
the eleven Archimedean tiling lattices (in which all vertices are of the same type), 
designated as (3,6,3,6) in the notation of [1], which means that each vertex touches a 
triangle, hexagon, triangle, and hexagon. The Kagome lattice is intimately related 
to other important lattices in percolation; its sites correspond to the bonds of the 
honeycomb lattice, which implies that the percolation thresholds Pc for those two 
systems are the same, 1 — 2 sin(7r/18), and by duality they are also equal to one minus 
the threshold for bond percolation on the triangular lattice [2] . 

For bond percolation on the Kagome lattice, no exact expression for is known, 
although two conjectures were made a number of years ago, both within the larger 
context of the g-state Potts model from which percolation follows in the limit of g = 1. 
Wu [3] conjectured that Pc is the solution to 

/ - 6/ + 12/ - - 3/ + 1 = (1) 

which yields 

p^(Wu) = 0.524429 718 . (2) 

Subsequently, Enting and Wu [4] showed that the general conjecture is not valid for 
g = 3, thereby casting doubt on its validity for g = 1. Tsallis [5] conjectured that Pc 
satisfies 

-p^ +J3 = 1 - 2sin(7r/18) (3) 

which yields 

j9c(Tsallis) = 0.522 372 078 . (4) 

Note that the quantity on the right-hand side of (3) is identical to the threshold 
for site percolation on the Kagome lattice and is the solution to the cubic equation 
— 3y^ -h 1 = 0, so it follows that (4) is the solution to the ninth order equation 

p^ - 3/ + d>p^ - 6p^ - + 5p^ + 3p^ - 1 = . (5) 

The only existing numerical values of Pc of relatively high precision appear to be 
those of Yonezawa et al. [6], who find 0.5244±0.0002, and van der Marck [7], who finds 
0.5243 ± 0.0004. These results clearly favor Wu's value over Tsallis'. Note that Hu, 
Chen and Wu [8] have also recently presented numerical evidence that Wu's conjecture 
still works quite well (and better than Tsallis') for the Potts model of various q. In 
order to investigate the validity of these conjectures further, and to provide an accurate 
value of Pc for use by others [9] , we have carried out a new numerical study to determine 
the percolation threshold for bond percolation on the Kagome lattice. 
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2. Method 

The method we employed is the hull-gradient method [10], in which the gradient- 
percolation frontier is created by a hull-generating walk. In gradient percolation 
[11,12], a linear gradient in p is imposed on the lattice in the vertical direction; as 
the height increases, the occupied bond density p also increases. The estimate of 
Pc is related to the average position of the frontier of the percolating region. To 
simultaneously create and measure that frontier, a hull-generating walk is employed 
[13,14]. In this walk, the status of a bond (whether occupied or vacant) is determined 
when the bond is visited by generating a random number and comparing that number 
to the occupation probability for that height. 

The efficiency of this method derives from the fact that in the hull-generating 
method, the entire lattice is not ffiled before the walk begins. Rather, the lattice is 
initialized with all bonds undetermined, and the state of the bonds are decided only 
when they are visited. If the walk does not reach a given bond, then the status of that 
bond remains undetermined, and no random number need be generated. The error of 
this method is at the statistical limit 0.5A^^^/^, where is the total quantity of random 
numbers generated. Previously, we used the hull-gradient method to find Pc for site 
percolation on the square lattice to six significant figures, 0.592 746 ± 0.000 000 5 
[10,15]. This value was confirmed using a different method (also implemented using 
hulls) [16], and is consistent with the most precise value 0.592 77±0.000 05 [17] obtained 
using the more traditional average crossing-probability method [18]. In [10], we also 
applied the hull-gradient method to determine pc for site percolation on the Kagome 
lattice, and found a value (0.652 704 ± 0.000 009) in agreement with the theoretical 
prediction mentioned above. 

For bond percolation, the hull walk simplifies to a trajectory that "bounces" back 
and forth between the centers of the occupied and vacant bonds of the hull, as first 
noted by Grassberger [19] for the case of bond percolation on the square lattice. 
For the Kagome lattice, a similar method can be used. As shown in Fig. 1, the walk 
moves along line segments that connect centers of adjacent bonds. These line segments 
produce a new lattice whose topology is the Archimedean (6,4,3,4)-lattice [1] shown in 
Fig. 2. Here, the walker turns clockwise when an occupied bond is hit, and counter- 
clockwise when a vacant bond is hit, so that the bonds are effectively rotators [20] or 
mirrors [21,22] on the vertices of the (6,4,3,4)-lattice. An occupied bond (probability 
p) on the Kagome lattice corresponds to a mirror placed tangent to the vertex of the 
hexagon on the (6,4,3,4) lattice, while a vacant bond (probability 1—p) corresponds to 
a mirror that intersects the hexagon. The hull walk is then a mirror-model trajectory 
[21,22] on the (6,4,3,4)-lattice. 

Many other representations of the hull walk can also be made. The vacant bonds 
on the Kagome lattice can be associated with occupied bonds on the dice (or "diced" ) 
lattice shown in Fig. 3, which is dual to the Kagome lattice, and the walk creates a hull 
on that lattice also. The hull trajectory is also a self-avoiding trail on the directed 
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(6,4,3,4)-lattice, with opposing direction vectors at each vertex. The hulls on this 
lattice can also be produced by a random tiling similar to [23], with "kite" -shaped 
tiles having weights p and 1 — p as shown in Fig. 4. 

In gradient percolation, the hull of the percolating region resides on bonds whose 
average value of p gives an estimate Pc(fi'), which approaches as the gradient g = \Vp\ 
goes to zero [11,12]. Equivalent to taking the average value of p on bonds of the hull, 
one can simply take as the estimate [12] 

Pcig) = ^ ^^^l^ — , (6) 
"-occ + ?^vac 

since the expected fraction of occupied bonds in the hull equals the average value 
of p of the vertices that belong to the hull. Here, nocc is the number of vertices 
corresponding to occupied bonds and nyac is the number of vertices corresponding to 
vacant bonds in the hull. 

To represent the (6,4,3,4)-lattice in the computer efficiently, it is necessary 
transform it to align on a rectilinear grid. One way to do this is shown in Fig. 5, 
where the basic rectangle of six sites is repeated in every second column and every 
third row. While we have distorted the lattice laterally to accommodate the square 
lattice periodicity, we have not shifted any vertices vertically with respect to each 
other, so as not to affect the gradient in the vertical direction. In producing that 
gradient, we made the change in p proportional to the actual height, so that changes 
between the wider rows in Fig. 2 equal twice the change of the narrower ones. The 
gradient g is defined here hj g = Ap/i where Ap is the change of p between the wider 
rows, and £ is the bond length, taken to be unity. (We also considered two alternate 
representations where the gradient was not precisely uniform on a local scale; the 
behavior of these systems is discussed in the Appendix.) A 3d array was used to store 
the six possible outgoing directions based upon the two incoming directions, the six 
types of vertices, and the status (occupied or vacant) of the vertex. 

The lattice was initialized by filling the first column halfway with occupied bonds 
and the rest with vacant bonds, which prevents the walk from closing on itself at the 
start. With the gradient in the vertical direction, the walk naturally drifts to the right. 
Periodic boundary conditions were applied in the horizontal direction, and each new 
column to the right was cleared off as it was first visited. This allows the simulation 
to run indefinitely and have essentially no boundary effects from the horizontal ends 
of the system. We tracked the maximum distance the walk wandered to the left of 
the moving front in order to confirm that the system width was sufficient to preclude 
wraparound errors. 

To rule out systematic errors related to random number generation, three different 
generators were tried. For most of the runs, we used the shift-register sequence 
generator R7(9689) [16, 24] defined by 

Xfi = Xn-47l ^a;„_i586 ^a;„_6988 ^a;„_9689 (7) 

where ^ is the bitwise exclusive-or operation. This "four-tap" generator is equivalent 
to decimating by 7 (taking every seventh term) of the sequence generated by the 
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two-tap rule R(9689), Xn — Xn^^n ^x^.gesg, which follows from [25]. This decimation 
has the effect of vastly reducing the three- and four-point correlations of the two-tap 
generator, and was previously found to yield good behavior for problems of this type 
[26]. For the system with the second smallest gradient, we considered two additional 
random number generators. The second generator, R2i(9689), was obtained by further 
decimating R7(9689) three times, simply by using every third number, which may yield 
better statistical properties. For the third generator, we used a traditional congruential 
generator CONG, but with a very large modulus of 64 bits [27]: 

Xn = (5 081 641 266 417 562 522a;„_i + 11) mod 2^^ . (8) 

First, to study the general finite-size behavior, we considered lattices of height 
H = 128, 256, 512, 1024, 2048 and 4096, with widths sufficient to avoid wraparound 
error, and p ranging all the way from to 1. For these systems, g = 1.5/ H, the factor 
of 1.5 resulting from the changing increment of p between the wide and narrow rows 
as described above. Approximately 10^^ occupied plus vacant vertices were generated 
for each of these lattices. Then, to obtain a precise final value, we used systems with 
very small gradients g = 2.564 ■ 10^"^ and 7.324 • 10~^ by using lattices of height 4096 
and 8192, with p ranging from 0.49 to 0.56 and 0.505 to 0.545, respectively; 2 • 10^^ 
steps were carried out for each of these systems. In all, several months of workstation 
computer time were used to obtain the final results. 

In the simulations, we kept track of the maximum and minimum heights of the 
hull. As g decreases, the relative width of the walk decreases as g^^'^ [12], allowing us 
to expand the gradient as mentioned above. Note that we also carried out runs on a 
system of height 64, but ran into difficulty because the walk wandered all the way to 
a top or bottom boundary where it got stuck in a dead end. Presumably, this could 
be averted by more carefully constructing those boundaries, but we did not attempt 
to do it. 



3. Results and Discussion 

First we compare the three random number generators. Table 1 gives the results 
for runs for g = 2.564 ■ 10~^, where all generators were used. Error bars represent 
one standard deviation, and follow from the statistical formula [Pc(l ~ Pc) /N]^^"^ ~ 
0.5A^~^/^ where N — riocc + n^ac, since the occupancy of each vertex is achieved with 
complete statistical independence. Clearly, the three random number generators give 
statistically consistent results, and we thus averaged their results to get the final data 
point for this g. 

Our complete results are shown in Fig. 6, where we plot pdg) vs. g for all the 
lattices we considered. This plot provides good evidence that the dependence of Pc{g) 
upon g is linear (as observed by Rosso et al. [12] for the case of site percolation on the 
square lattice) with the behavior 

Pc{g)=Pc + 0.010 g (9) 
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The slope 0.010 implies that the average value of p differs from by an amount 
corresponding to one hundredth of a lattice spacing, independent of the gradient. For 
g = 2.564 • 10~^, this implies a finite-size correction Pc{g) — Pc = 2.3 • 10~^ that is 
somewhat smaller than the error bars given in Table 1. For g = 7.324 • 10~^, the 
simulations of 2.2 • 10^^ steps using R7(9689) yielded 0.524405 5 ±3.4 ■ 10"^. Here, the 
finite-size correction is insignificant compared to the statistical error. Putting these 
results all together, we obtain our final result 

Pe = 0.524 405 3 ±0.000 000 3 (10) 

This result is consistent with — but nearly 1000 times more precise than — previous 
values [6,7]. It evidently agrees with neither Tsallis' nor Wu's predictions, although 
the difference with Wu's approximate conjecture is remarkably small, only 47 parts 
per million (but still much larger than our error bars of less than one part per million) . 
Note that it would be quite difficult to observe this small difference in pc using 
conventional methods [e.g., 6, 7, 8, 17, 18, 28, 29]. 

If neither Wu's nor Tsallis' conjecture is valid, is there perhaps some other simple 
polynomial that yields Pc? In the absence of a theory, we can search numerically for 
possible candidates consistent with our numerical result. However, if we allow the 
maximum order of the polynomial to be six, and the integer coefficients to be as large 
as say ±24 (except for the leading coefficient, which we restrict to unity), then we find 
literally thousands of polynomials with roots within two standard deviations of (10). 



Some examples are: 

p^±7p^±17p-10 = , = 0.524405 335 (11a) 

/-24p2±p±6 = , pc = 0.524405 671 (lib) 

p^-2j9^-2p3±16p^-4 = , Pc = 0.524405 424 (11c) 

p^±5p=^-8p2±18p-8 = , pc = 0.524404 863 (lid) 

p^ + 3p^ - 3p^ ± 12p - 6 = , Pc = 0.524405 290 (lie) 

p^ ± 5p^ ± 7p^ - 5p2 ± 6p - 3 = , Pc = 0.524 405 306 (llf ) 

p^ ± 3p^ ± 9p^ - p^ ± p^ ± 2p - 2 = , Pc = 0.524 405 134 (llg) 



Note also that (11/40)^/^ = 0.524 404 424 is only slightly low. Unfortunately, it is not 
possible by numerical means to determine Pc with sufficient accuracy to distinguish 
which of these many polynomials is the correct one (if indeed one is!). 

Note finally that (10) implies that the bond threshold of the dice lattice shown in 
Fig. 3 is given by 

Pc(dice) = 1 -pc(Kagome) = 0.475 594 7 ± 0.000 0003 (12) 
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Table 1. Results for p^{g) given by (6) for runs with height H = 4096 and gradient 
g = 0.000 025 64, using three different random number generators (RNG). is the 
total number of occupied and vacant bonds generated, and a represents one standard 
deviation of error (68% confidence interval). 



RNG 


N 




pM 


£7 = 0.5iV-°-5 


R7(9689) 


1.0- 


10^2 


0.5244048 


±0.000000 5 


R2i(9689) 


0.5- 


10^2 


0.524405 3 


±0.000000 7 


CONG (64-bit) 


0.5 • 


1012 


0.524 405 9 


±0.000 000 7 


Average 


2.0 ■ 


1012 


0.524405 2 


±0.000 0004 



8 Robert M. Ziff and Paul N. Suding 

Appendix. 

Besides the system described above with the gradient apphed completely 
uniformly, we also considered two systems in which gradient was not constructed 
so precisely on a local scale, and it is instructive to see their effects on the finite-size 
behavior. At first, we squared-off the (6,3,4,3)-lattice by "stretching" it horizontally, 
leading to a lattice similar to Fig. 5 but rotated by 90°. Thus, we effectively pushed 
up and down alternating columns in the original lattice, and the gradient was applied 
equally between all the rows in this distorted lattice. The idea was that these local 
variations should have little effect on the behavior when the gradient is small. However, 
the deviations turned out to be rather large, until the gradient dropped to about 0.002, 
as shown in Fig. 7 (case A). As a second test (case B in Fig. 7), we represented the 
lattice as in Fig. 5, but applied the gradient equally between all rows, whether "wide" 
or "narrow." Again, the behavior of the finite-size corrections to Pc(fi') differed from 
the uniform case, but not as much here. In the limit of g small, where p changes 
very little from row to row, all systems followed the same limiting finite-size behavior 
as (9). These results show, however, that for the linear behavior to remain valid for 
moderately large gradients, a uniform gradient must be applied to the lattice in its 
actual configuration. 

Acknowledgments 

The authors thank A J Guttmann for interesting them in this problem, and thank 
F. Y. Wu for comments. This material is based upon work supported by the U. S. 
National Science Foundation under Grant No. DMR-9520700. 

References 

E-mail addresses: rziff@umich.edu, psuding@umich.edu 

[1] Griinbaum B and Shephard G C 1987 Tilings and Patterns (New York: Freeman and Co.) 

[2] Sykcs M F and Essam J W 1964 J. Math. Phys. 5 1117 

[3] Wu F Y 1979 J. Phys. C: Solid State 12 L645 

[4] Enting I G and Wu F Y 1982 J. Stat. Phys. 28 351 

[5] TsaUis C 1982 J. Phys. C: Solid State 15 L757 

[6] Yonezawa F, Sakamoto S and Hori M 1989 Phys. Rev. B 40 636 

[7] van der Marck S C 1997 Phys. Rev. E 55 1514 

[8] Hti C-K, Chen J-A, and Wu F Y 1994 Mod. Phys. Lett. B 8 455 

[9] Guttmann A J and Jensen I 1997 (to be pubhshed) 

[10] Ziff R M and Sapoval B 1986 J. Phys. A: Math. Gen. 18 L1169 

[11] Sapoval B, Rosso M and Gouyct J F 1985 J. Physique Lett. (Paris) 46 L149 

[12] Rosso M, Gouyet J F and Sapoval B 1986 Phys. Rev. B 32 6053 

[13] Ziff R M, Cummings P T and Stell G 1984 J. Phys. A: Math. Gen. 17 3009 

[14] Wcinrib A and Trugman S 1985 Phys. Rev. B 31 2993 

[15] Ziff R M and Stell G 1988 UM LaSC Report 88-4 (unpublished) 

[16] Ziff R M 1992 Phys. Rev. Lett. 69 2670 

[17] Gebele T 1984 J. Phys. A: Math. Gen. 17 L51 



Bond percolation threshold on the Kagome lattice 



9 



[18] Stauffer D and Aharony A 1994 An Introduction to Percolation Theory, 2nd revised edition 

(London: Taylor and Francis) 

[19] Grassberger P 1984 J. Phys. A: Math. Gen. 19 2675 

[20] Gunn J M F and Ortuno M 1985 J. Phys. A: Math. Gen. 18 L1096 

[21] Cao M-S and Cohen E G D 1997 J. Statis. Phys. (to be published) 

[22] Wang F and Cohen E G D 1995 J. Statis. Phys. 81 467 

[23] Roux S, Guyon E and Sornette D 1988 J. Phys. A: Math. Gen. 21 L475 

[24] Golomb S W 1967 Shift Register Sequences (San Francisco: Holden-Day) 

[25] Zierler N 1969 Inform. Control 15 67 

[26] ZiffRM unpublished 

[27] Wood W W (Los Alamos, NM) private communication 

[28] Reynolds P J, Stanley H E and Klein W 1980 Phys. Rev. B 21 1223 

[29] Hu C-K 1992 Phys. Rev. Lett. 69 2739 



10 



Robert M. Ziff and Paul N. Suding 



Figure captions 

Figure 1. The Kagome lattice, showing a path of bonds that represent the occupied 
bonds of the frontier of the percolating region, which is above it (gradient increases 

in the vertical direction). The dashed line shows the hull walk, which "bounces" back 
and forth between centers of the occupied and vacant bonds of the hull. 

Figure 2. The (6,4,3,4)-Archimcdcan lattice, on which the hull walk of Fig. 2 effectively 
takes place. The trajectory of the walk is equivalent to a mirror-model trajectory [21,22] 
on the (6,4,3,4)-lattice, in which the bonds on the underlying Kagome lattice act as 
the mirrors. 

Figure 3. The same hull as in Fig. 2 on the dice (or "diced") lattice, the dual to the 

Kagome lattice. Vacant bonds on the Kagome lattice correspond to occupied bonds on 
the dice, and p^(dice) = 1 —p^, (Kagome). 

Figure 4. A tiling representation of the hull walk on the Kagome lattice, shown in the 
central part of the figure. The tile marked with probability p corresponds to an occupied 
bond on the underlying Kagome lattice, while the one marked 1—p corresponds to a 
vacant bond. 

Figure 5. The representation of the (6,4,3,4)-lattice on a rectilinear grid, as utilized 
in the computer program's two-dimensional array. 

Figure 6. A plot of the estimate Pc{g) determined by (6), versus the gradient g, 
implying the linear relation (9). The error bars represent one standard deviation of 
statistical error. 

Figure 7. A similar plot as in Fig. 6, with data from two systems (case A, dashed line 
through data points, and case B, solid line through data points) in which the gradient 
is not locally uniform, as described in the Appendix. The straight line represents the 
fit of the data of Fig. 6 for the system the uniform gradient. For small g, all systems 
follow the same asymptotic behavior. 



